! Bessel J_1(x) function in double precision
!
      function dbesj1(x)
      implicit real*8 (a - h, o - z)
      dimension a(0 : 7), b(0 : 64), c(0 : 69), d(0 : 51)
      parameter (pi4 = 0.78539816339744830962d0)
      data (a(i), i = 0, 7) / 
     &    -0.00000000000014810349d0, 0.00000000003363594618d0, 
     &    -0.00000000565140051697d0, 0.00000067816840144764d0, 
     &    -0.00005425347222188379d0, 0.00260416666666662438d0, 
     &    -0.06249999999999999799d0, 0.49999999999999999998d0 / 
      data (b(i), i = 0, 12) / 
     &    0.00000000000243721316d0, -0.00000000009400554763d0, 
     &    0.00000000306053389980d0, -0.00000008287270492518d0, 
     &    0.00000183020515991344d0, -0.00003219783841164382d0, 
     &    0.00043795830161515318d0, -0.00442952351530868999d0, 
     &    0.03157908273375945955d0, -0.14682160488052520107d0, 
     &    0.39309619054093640008d0, -0.47952808215101070280d0, 
     &    0.14148999344027125140d0 / 
      data (b(i), i = 13, 25) / 
     &    0.00000000000182119257d0, -0.00000000006862117678d0, 
     &    0.00000000217327908360d0, -0.00000005693592917820d0, 
     &    0.00000120771046483277d0, -0.00002020151799736374d0, 
     &    0.00025745933218048448d0, -0.00238514907946126334d0, 
     &    0.01499220060892984289d0, -0.05707238494868888345d0, 
     &    0.10375225210588234727d0, -0.02721551202427354117d0, 
     &    -0.06420643306727498985d0 / 
      data (b(i), i = 26, 38) / 
     &    0.000000000001352611196d0, -0.000000000049706947875d0, 
     &    0.000000001527944986332d0, -0.000000038602878823401d0, 
     &    0.000000782618036237845d0, -0.000012349994748451100d0, 
     &    0.000145508295194426686d0, -0.001203649737425854162d0, 
     &    0.006299092495799005109d0, -0.016449840761170764763d0, 
     &    0.002106328565019748701d0, 0.058527410006860734650d0, 
     &    -0.031896615709705053191d0 / 
      data (b(i), i = 39, 51) / 
     &    0.000000000000997982124d0, -0.000000000035702556073d0, 
     &    0.000000001062332772617d0, -0.000000025779624221725d0, 
     &    0.000000496382962683556d0, -0.000007310776625173004d0, 
     &    0.000078028107569541842d0, -0.000550624088538081113d0, 
     &    0.002081442840335570371d0, -0.000771292652260286633d0, 
     &    -0.019541271866742634199d0, 0.033361194224480445382d0, 
     &    0.017516628654559387164d0 / 
      data (b(i), i = 52, 64) / 
     &    0.000000000000731050661d0, -0.000000000025404499912d0, 
     &    0.000000000729360079088d0, -0.000000016915375004937d0, 
     &    0.000000306748319652546d0, -0.000004151324014331739d0, 
     &    0.000038793392054271497d0, -0.000211180556924525773d0, 
     &    0.000274577195102593786d0, 0.003378676555289966782d0, 
     &    -0.013842821799754920148d0, -0.002041834048574905921d0, 
     &    0.032167266073736023299d0 / 
      data (c(i), i = 0, 13) / 
     &    -0.00000000001185964494d0, 0.00000000039110295657d0, 
     &    0.00000000180385519493d0, -0.00000005575391345723d0, 
     &    -0.00000018635897017174d0, 0.00000542738239401869d0, 
     &    0.00001181490114244279d0, -0.00033000319398521070d0, 
     &    -0.00037717832892725053d0, 0.01070685852970608288d0, 
     &    0.00356629346707622489d0, -0.13524776185998074716d0, 
     &    0.00980725611657523952d0, 0.27312196367405374425d0 / 
      data (c(i), i = 14, 27) / 
     &    -0.00000000003029591097d0, 0.00000000009259293559d0, 
     &    0.00000000496321971223d0, -0.00000001518137078639d0, 
     &    -0.00000057045127595547d0, 0.00000171237271302072d0, 
     &    0.00004271400348035384d0, -0.00012152454198713258d0, 
     &    -0.00184155714921474963d0, 0.00462994691003219055d0, 
     &    0.03671737063840232452d0, -0.06863857568599167175d0, 
     &    -0.21090395092505707655d0, 0.16126443075752985095d0 / 
      data (c(i), i = 28, 41) / 
     &    -0.00000000002197602080d0, -0.00000000027659100729d0, 
     &    0.00000000374295124827d0, 0.00000003684765777023d0, 
     &    -0.00000045072801091574d0, -0.00000327941630669276d0, 
     &    0.00003571371554516300d0, 0.00017664005411843533d0, 
     &    -0.00165119297594774104d0, -0.00485925381792986774d0, 
     &    0.03593306985381680131d0, 0.04997877588191962563d0, 
     &    -0.22913866929783936544d0, -0.07885001422733148814d0 / 
      data (c(i), i = 42, 55) / 
     &    0.00000000000516292316d0, -0.00000000039445956763d0, 
     &    -0.00000000066220021263d0, 0.00000005511286218639d0, 
     &    0.00000005012579400780d0, -0.00000522111059203425d0, 
     &    -0.00000134311394455105d0, 0.00030612891890766805d0, 
     &    -0.00007103391195326182d0, -0.00949316714311443491d0, 
     &    0.00455036998246516948d0, 0.11540391585989614784d0, 
     &    -0.04779493761902840455d0, -0.22837862066532347460d0 / 
      data (c(i), i = 56, 69) / 
     &    0.00000000002697817493d0, -0.00000000016633326949d0, 
     &    -0.00000000433134860350d0, 0.00000002508404686362d0, 
     &    0.00000048528284780984d0, -0.00000258267851112118d0, 
     &    -0.00003521049080466759d0, 0.00016566324273339952d0, 
     &    0.00146474737522491617d0, -0.00565140892697147306d0, 
     &    -0.02833882055679300400d0, 0.07580744376982855057d0, 
     &    0.16012275906960187978d0, -0.16548380461475971845d0 / 
      data (d(i), i = 0, 12) / 
     &    -1.272346002224188092d-14, 3.370464692346669075d-13, 
     &    -1.144940314335484869d-11, 6.863141561083429745d-10, 
     &    -9.491933932960924159d-8, 5.301676561445687562d-5, 
     &    0.1628675039676399740d0, -3.652982212914147794d-13, 
     &    1.151126750560028914d-11, -5.165585095674343486d-10, 
     &    4.657991250060549892d-8, -1.186794704692706504d-5, 
     &    1.562499999999994026d-2 / 
      data (d(i), i = 13, 25) / 
     &    -8.713069680903981555d-15, 3.140780373478474935d-13, 
     &    -1.139089186076256597d-11, 6.862299023338785566d-10, 
     &    -9.491926788274594674d-8, 5.301676558106268323d-5, 
     &    0.1628675039676466220d0, -2.792555727162752006d-13, 
     &    1.108650207651756807d-11, -5.156745588549830981d-10, 
     &    4.657894859077370979d-8, -1.186794650130550256d-5, 
     &    1.562499999987299901d-2 / 
      data (d(i), i = 26, 38) / 
     &    -6.304859171204770696d-15, 2.857249044208791652d-13, 
     &    -1.124956921556753188d-11, 6.858482894906716661d-10, 
     &    -9.491867953516898460d-8, 5.301676509057781574d-5, 
     &    0.1628675039678191167d0, -2.185193490132496053d-13, 
     &    1.048820673697426074d-11, -5.132819367467680132d-10, 
     &    4.657409437372994220d-8, -1.186794150862988921d-5, 
     &    1.562499999779270706d-2 / 
      data (d(i), i = 39, 51) / 
     &    -4.740417209792009850d-15, 2.578715253644144182d-13, 
     &    -1.104148898414138857d-11, 6.850134201626289183d-10, 
     &    -9.491678234174919640d-8, 5.301676277588728159d-5, 
     &    0.1628675039690033136d0, -1.755122057493842290d-13, 
     &    9.848723331445182397d-12, -5.094535425482245697d-10, 
     &    4.656255982268609304d-8, -1.186792402114394891d-5, 
     &    1.562499998712198636d-2 / 
      w = abs(x)
      if (w .lt. 1) then
          t = w * w
          y = (((((((a(0) * t + a(1)) * t + 
     &        a(2)) * t + a(3)) * t + a(4)) * t + 
     &        a(5)) * t + a(6)) * t + a(7)) * w
      else if (w .lt. 8.5d0) then
          t = w * w * 0.0625d0
          k = int(t)
          t = t - (k + 0.5d0)
          k = k * 13
          y = ((((((((((((b(k) * t + b(k + 1)) * t + 
     &        b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + 
     &        b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + 
     &        b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + 
     &        b(k + 11)) * t + b(k + 12)) * w
      else if (w .lt. 12.5d0) then
          k = int(w)
          t = w - (k + 0.5d0)
          k = 14 * (k - 8)
          y = ((((((((((((c(k) * t + c(k + 1)) * t + 
     &        c(k + 2)) * t + c(k + 3)) * t + c(k + 4)) * t + 
     &        c(k + 5)) * t + c(k + 6)) * t + c(k + 7)) * t + 
     &        c(k + 8)) * t + c(k + 9)) * t + c(k + 10)) * t + 
     &        c(k + 11)) * t + c(k + 12)) * t + c(k + 13)
      else
          v = 24 / w
          t = v * v
          k = 13 * (int(t))
          y = ((((((d(k) * t + d(k + 1)) * t + 
     &        d(k + 2)) * t + d(k + 3)) * t + d(k + 4)) * t + 
     &        d(k + 5)) * t + d(k + 6)) * sqrt(v)
          theta = (((((d(k + 7) * t + d(k + 8)) * t + 
     &        d(k + 9)) * t + d(k + 10)) * t + d(k + 11)) * t + 
     &        d(k + 12)) * v - pi4
          y = y * sin(w + theta)
      end if
      if (x .lt. 0) y = -y
      dbesj1 = y
      end
!
